function [v,inten,pi,p,q,sigma,alpha] = HJB_diffusion_implicit(A,K,...
    alpha,sigma,gamma,rho,maxit,crit,I,J)

    % state space
    grid_p = linspace(0,1,I);
    dp = 1/(I-1);

    grid_q = linspace(0,1,J);
    dq = 1/(J-1);

    states = combvec(grid_p,grid_q)';
    p = states(:,1); 
    q = states(:,2);

    % investment value
    pi = max(0, -K + A*2*alpha*(p-q)); % investment value
    v0 = pi;
    
    % Placeholders
    p = reshape(p,I,J);
    q = reshape(q,I,J);
    v0 = reshape(v0,I,J);
    pi = reshape(pi,I,J);
    
    xi = 2*alpha/sigma;

    % Value function iteration (given i,j)
    v = v0;
    
    for n=1:maxit
        V = v;

        % Computation of optimal intensity given V
        inten = intensity(p,dp,q,dq,V,alpha,sigma,gamma);

        % Coef in front of deriv terms
        Sq = q.*(1-q).*xi.*(2*alpha*(p-q)./sigma);
        Sqq = (q.^2).*((1-q).^2)*xi^2;    
        Spp = ((p.^2).*((1-p).^2)*xi^2).*(1+inten); 
        Spq = p.*(1-p).*q.*(1-q).*sqrt(1+inten)*xi^2;

        % Vectorized HJB equation coefficients 
        Xhat = -min(Sq,0)/dq + Sqq/(2*dq^2); % V_{i,j-1} coef
        Yhat = min(Sq,0)/dq - max(Sq,0)/dq - Sqq/dq^2; % V_{i,j} coef
        Zhat = max(Sq,0)/dq + Sqq/(2*dq^2); % V_{i,j+1} coef
        X = Spp/(2*dp^2); % V_{i-1,j} coef
        Y = -Spp/dp^2; % V_{i,j} coef
        Z = Spp/(2*dp^2); % V_{i+1,j} coef
        Xbar = Spq/(8*dp*dq); % V_{i-1,j-1} coef
        Ybar = -Spq/(8*dp*dq); % V_{i-1,j+1} coef
        Zbar = -Spq/(8*dp*dq); % V_{i+1,j-1} coef
        Kbar = Spq/(8*dp*dq); % V_{i+1,j+1} coef

        % Adjustments for Reflecting Boundaries
        ijndiag = Xhat(:); ijdiag = Y(:) + Yhat(:); ijpdiag = Zhat(:);
        injndiag = Xbar(:); injdiag = X(:); injpdiag = Ybar(:);
        ipjndiag = Zbar(:); ipjdiag = Z(:); ipjpdiag = Kbar(:);

        % Sparse HJB representation, rV = -c + AV
        A = sparse_matrix_function(injndiag,ijndiag,ipjndiag,injdiag,ijdiag,...
            ipjdiag,injpdiag,ijpdiag,ipjpdiag,I,J);      

        % Solve problem as LCP
        B = rho*speye(I*J) - A;
        m = inten(:).^gamma + B*pi(:); 
        z0 = V(:)-pi(:);

        l = zeros(I*J,1); u = max(max(pi))*ones(I*J,1);
        z = LCP(B,m,l,u,z0,0);

        V_stack = z+pi(:);    
        v = reshape(V_stack,I,J);

        Vchange = v - V;
        dist = max(max(abs(Vchange)));

        out = ['Iteration ',num2str(n),', Max Difference ',num2str(dist)];
        disp(out) % iteration number
        if dist<crit
            out = ['Value Function Converged at Iteration ',num2str(n)];
            disp(out)
            break
        end
    end